Run 5.6 PHOTOSYNTHESIS Analysis, Script Chunks, and Plots

This is the analysis of the final run of the Ulva lactuca salinity and nutrient experiments conducted on the lanai in St. John 616. These experiments incorporated four paired salinity and nutrient levels with three temperature levels. Each run produced an n = 2 and was repeated initially 8 times for a total of n = 16. Data gaps were identified and filled by early March 2022. This output reflects all data totaling five treatments for Ulva lactuca.

Packages loaded:

library(lme4)
library(lmerTest)
library(effects)
library(car)
library(MuMIn)
library (dplyr)
library(emmeans)
library(DHARMa)
library(performance)
library(patchwork)
library(ggplot2)
library(ggpubr)

Load and prepare the dataset

Open the output dataset generated by the ps_script_clean_to_ek_alpha.R script in the phytotools_alpha_ek project

run5_6_photosyn_data <- read.csv("/Users/Angela/src/work/limu/algal_growth_photosynthesis/data_input/run5-6_ek_alpha.csv")

Assign run as a factor

run5_6_photosyn_data$Run <- as.factor(run5_6_photosyn_data$Run)

Assign temperature as a factor

run5_6_photosyn_data$Temperature <- as.factor(run5_6_photosyn_data$Temp...C.)

Assign treatment as characters from integers then to factors

run5_6_photosyn_data$Treatment <- as.factor(as.character(run5_6_photosyn_data$Treatment))

Subset the data and toggle between the species for output. Use Day 9 for final analysis ONLY

ulva <- subset(run5_6_photosyn_data, Species == "ul" & RLC.Day == 9)

Add a column for growth rate from growth rate dataset to the alredy subsetted ulva data frame

growth_rate <- read.csv("/Users/Angela/src/work/limu/algal_growth_photosynthesis/data_input/run5-6_growth_all_042922.csv")
gr_ulva <- subset(growth_rate, Species == "Ul")
ulva$growth_rate <- round((gr_ulva$final.weight - gr_ulva$Initial.weight) / gr_ulva$Initial.weight * 100, digits = 2)

rETRmax – Run the model

Run model for rETRmax with two fixed effect variables and three random effects variables

run5_6_photosyn_model_noint <- lmer(formula = rETRmax ~ Treatment + Temperature 
                                    + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), 
                                    data = ulva)

rETRmax – Make a histogram and residual plots of the data

hist(ulva$rETRmax, main = paste("Ulva lactuca rETRmax"), col = "darkolivegreen2", labels = TRUE)

plot(resid(run5_6_photosyn_model_noint) ~ fitted(run5_6_photosyn_model_noint))

qqnorm(resid(run5_6_photosyn_model_noint))
qqline(resid(run5_6_photosyn_model_noint))

rETRmax – Check the performance of the model

performance::check_model(run5_6_photosyn_model_noint)

These outputs show the model is acceptable

rETRmax – Check r2 for model fit and print the model statistics summary

r.squaredGLMM(run5_6_photosyn_model_noint)
## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.
##            R2m       R2c
## [1,] 0.5882889 0.6835032
summary(run5_6_photosyn_model_noint)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: rETRmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) +  
##     (1 | RLC.Order)
##    Data: ulva
## 
## REML criterion at convergence: 1827.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8253 -0.6047  0.0845  0.5200  3.4082 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Plant.ID  (Intercept)  22.155   4.707  
##  RLC.Order (Intercept)   5.044   2.246  
##  Run       (Intercept)   6.305   2.511  
##  Residual              111.367  10.553  
## Number of obs: 240, groups:  Plant.ID, 95; RLC.Order, 6; Run, 5
## 
## Fixed effects:
##               Estimate Std. Error     df t value Pr(>|t|)    
## (Intercept)     34.819      3.437  3.609  10.130 0.000887 ***
## Treatment1      20.641      3.669  3.083   5.626 0.010318 *  
## Treatment2      21.727      3.669  3.083   5.922 0.008906 ** 
## Treatment3      37.848      3.669  3.083  10.317 0.001724 ** 
## Treatment4      39.172      3.669  3.083  10.677 0.001555 ** 
## Temperature27   -4.430      2.504 30.316  -1.770 0.086835 .  
## Temperature30   -2.157      2.293 65.525  -0.941 0.350296    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trtmn1 Trtmn2 Trtmn3 Trtmn4 Tmpr27
## Treatment1  -0.721                                   
## Treatment2  -0.721  0.828                            
## Treatment3  -0.721  0.828  0.828                     
## Treatment4  -0.721  0.828  0.828  0.828              
## Temperatr27 -0.348  0.003  0.003  0.003  0.003       
## Temperatr30 -0.338  0.000  0.000  0.000  0.000  0.475

rETRmax – Run an ANOVA and pairwise comparisons based on ANOVA results Results for temperature are non-significant, thus no pairwise comparisons for this

anova(run5_6_photosyn_model_noint, type = c("III"), ddf = "Satterthwaite")
ulva_photosyn_model_aov <- aov(rETRmax ~ Treatment + Temperature, data = ulva)
TukeyHSD(ulva_photosyn_model_aov, "Treatment", ordered = FALSE)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = rETRmax ~ Treatment + Temperature, data = ulva)
## 
## $Treatment
##          diff       lwr       upr     p adj
## 1-0 20.610833 13.932211 27.289455 0.0000000
## 2-0 21.696250 15.017628 28.374872 0.0000000
## 3-0 37.818125 31.139503 44.496747 0.0000000
## 4-0 39.141667 32.463045 45.820289 0.0000000
## 2-1  1.085417 -5.593205  7.764039 0.9917098
## 3-1 17.207292 10.528670 23.885914 0.0000000
## 4-1 18.530833 11.852211 25.209455 0.0000000
## 3-2 16.121875  9.443253 22.800497 0.0000000
## 4-2 17.445417 10.766795 24.124039 0.0000000
## 4-3  1.323542 -5.355080  8.002164 0.9824820

rETRmax – Plot the results

plot(allEffects(run5_6_photosyn_model_noint))

Plot a regression between the photosynthetic independent variables of interest and growth rate

ulva_growth_etr_graph <- ggplot(ulva, aes(x=rETRmax, y=growth_rate)) + geom_point() + 
  geom_smooth(method = "lm", col = "black") + theme_bw() + 
  labs(title = "Ulva lactuca rETRmax vs Growth Rate", x = "rETRmax (μmols electrons m-2 s-1)", 
       y = "growth rate (%)") + stat_regline_equation(label.x = 25, label.y = 165) + stat_cor()
ulva_growth_etr_graph
## `geom_smooth()` using formula 'y ~ x'

Temperature did not have a significant effect on the outcome for rETRmax but it did for Ek and alpha. alpha histogram is NOT Normal, probably not a good fit for this analysis

Ek – Run the model

Run model for minimum saturating irradiance (Ek) with two fixed effect variables and three random effects variables

run5_6_photosyn_model_noint <- lmer(formula = ek.1 ~ Treatment + Temperature 
                                    + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), 
                                    data = ulva)

Ek – Make a histogram and residual plots of the data for ulva

hist(ulva$ek.1, main = paste("Ulva lactuca Ek"), col = "darkolivegreen3", labels = TRUE)

plot(resid(run5_6_photosyn_model_noint) ~ fitted(run5_6_photosyn_model_noint))

qqnorm(resid(run5_6_photosyn_model_noint))
qqline(resid(run5_6_photosyn_model_noint))

Ek – Check the performance of the model

performance::check_model(run5_6_photosyn_model_noint)

Ek – Check r2 for model fit and print the model statistics summary

r.squaredGLMM(run5_6_photosyn_model_noint)
##            R2m       R2c
## [1,] 0.5867698 0.7012116
summary(run5_6_photosyn_model_noint)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ek.1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) +  
##     (1 | RLC.Order)
##    Data: ulva
## 
## REML criterion at convergence: 2102.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.04017 -0.53517 -0.01583  0.51382  2.98662 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Plant.ID  (Intercept)  66.38    8.147  
##  RLC.Order (Intercept)  10.92    3.304  
##  Run       (Intercept)  61.94    7.870  
##  Residual              363.52   19.066  
## Number of obs: 240, groups:  Plant.ID, 95; RLC.Order, 6; Run, 5
## 
## Fixed effects:
##               Estimate Std. Error     df t value Pr(>|t|)   
## (Intercept)     16.796      8.856  3.225   1.897  0.14770   
## Treatment1      30.241      9.766  3.068   3.096  0.05187 . 
## Treatment2      40.947      9.766  3.068   4.193  0.02368 * 
## Treatment3      68.624      9.766  3.068   7.027  0.00550 **
## Treatment4      71.137      9.766  3.068   7.284  0.00495 **
## Temperature27   -9.826      4.311 29.565  -2.279  0.03004 * 
## Temperature30   -8.349      4.020 68.996  -2.077  0.04154 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trtmn1 Trtmn2 Trtmn3 Trtmn4 Tmpr27
## Treatment1  -0.820                                   
## Treatment2  -0.820  0.921                            
## Treatment3  -0.820  0.921  0.921                     
## Treatment4  -0.820  0.921  0.921  0.921              
## Temperatr27 -0.235  0.002  0.002  0.002  0.002       
## Temperatr30 -0.229  0.000  0.000  0.000  0.000  0.480

Ek – Run an ANOVA and pairwise comparisons based on ANOVA results

anova(run5_6_photosyn_model_noint, type = c("III"), ddf = "Satterthwaite")
ulva_photosyn_model_aov <- aov(ek.1 ~ Treatment + Temperature, data = ulva)
TukeyHSD(ulva_photosyn_model_aov, "Treatment", ordered = FALSE)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = ek.1 ~ Treatment + Temperature, data = ulva)
## 
## $Treatment
##         diff       lwr      upr     p adj
## 1-0 30.19583 17.895251 42.49642 0.0000000
## 2-0 40.90208 28.601501 53.20267 0.0000000
## 3-0 68.57917 56.278584 80.87975 0.0000000
## 4-0 71.09167 58.791084 83.39225 0.0000000
## 2-1 10.70625 -1.594333 23.00683 0.1208714
## 3-1 38.38333 26.082751 50.68392 0.0000000
## 4-1 40.89583 28.595251 53.19642 0.0000000
## 3-2 27.67708 15.376501 39.97767 0.0000000
## 4-2 30.18958 17.889001 42.49017 0.0000000
## 4-3  2.51250 -9.788083 14.81308 0.9803988
TukeyHSD(ulva_photosyn_model_aov, "Temperature", ordered = FALSE)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = ek.1 ~ Treatment + Temperature, data = ulva)
## 
## $Temperature
##          diff       lwr        upr     p adj
## 27-20 -9.9375 -18.11226 -1.7627403 0.0125118
## 30-20 -9.0075 -17.18226 -0.8327403 0.0267837
## 30-27  0.9300  -7.24476  9.1047597 0.9610888

Ek – Plot the results

plot(allEffects(run5_6_photosyn_model_noint))

ulva_growth_ek_graph <- ggplot(ulva, aes(x=ek.1, y=growth_rate)) + geom_point() + 
  geom_smooth(method = "lm", col = "black") + theme_bw() + 
  labs(title = "Ulva lactuca Ek vs Growth Rate", x = "Ek (μmols photons m-2 s-1)", y = "growth rate (%)") + 
  stat_regline_equation() + stat_cor(label.y = 160)
ulva_growth_ek_graph
## `geom_smooth()` using formula 'y ~ x'

alpha – Run the model

Run model for alpha withtwo fixed effect variables and three random effects variables

run5_6_photosyn_model_noint <- lmer(formula = alpha.1 ~ Treatment + Temperature 
                                    + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), 
                                    data = ulva)
## boundary (singular) fit: see help('isSingular')

alpha – Make a histogram and residual plots of the data for ulva

hist(ulva$alpha.1, main = paste("Ulva lactuca alpha"), col = "darkolivegreen4", labels = TRUE)

plot(resid(run5_6_photosyn_model_noint) ~ fitted(run5_6_photosyn_model_noint))

qqnorm(resid(run5_6_photosyn_model_noint))
qqline(resid(run5_6_photosyn_model_noint))

alpha – Check the performance of the model

performance::check_model(run5_6_photosyn_model_noint)

alpha – Check r2 for model fit and print the model statistics summary

r.squaredGLMM(run5_6_photosyn_model_noint)
##            R2m       R2c
## [1,] 0.4524773 0.5264424
summary(run5_6_photosyn_model_noint)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: alpha.1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) +  
##     (1 | RLC.Order)
##    Data: ulva
## 
## REML criterion at convergence: 453
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5531 -0.5754 -0.1204  0.3199  3.6187 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Plant.ID  (Intercept) 0.01181  0.1087  
##  RLC.Order (Intercept) 0.00000  0.0000  
##  Run       (Intercept) 0.04206  0.2051  
##  Residual              0.34486  0.5872  
## Number of obs: 240, groups:  Plant.ID, 95; RLC.Order, 6; Run, 5
## 
## Fixed effects:
##               Estimate Std. Error       df t value Pr(>|t|)   
## (Intercept)    2.15229    0.22956  3.28704   9.376  0.00175 **
## Treatment1    -0.84972    0.25969  3.44592  -3.272  0.03822 * 
## Treatment2    -1.11836    0.25969  3.44592  -4.306  0.01731 * 
## Treatment3    -1.48303    0.25969  3.44592  -5.711  0.00723 **
## Treatment4    -1.52872    0.25969  3.44592  -5.887  0.00656 **
## Temperature27  0.31167    0.09830 42.27696   3.171  0.00283 **
## Temperature30  0.25233    0.09804 46.29620   2.574  0.01332 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trtmn1 Trtmn2 Trtmn3 Trtmn4 Tmpr27
## Treatment1  -0.830                                   
## Treatment2  -0.830  0.893                            
## Treatment3  -0.830  0.893  0.893                     
## Treatment4  -0.830  0.893  0.893  0.893              
## Temperatr27 -0.214  0.001  0.001  0.001  0.001       
## Temperatr30 -0.214  0.000  0.000  0.000  0.000  0.499
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

alpha – Run an ANOVA and pairwise comparisons based on ANOVA results

anova(run5_6_photosyn_model_noint, type = c("III"), ddf = "Satterthwaite")
ulva_photosyn_model_aov <- aov(alpha.1 ~ Treatment + Temperature, data = ulva)
TukeyHSD(ulva_photosyn_model_aov, "Treatment", ordered = FALSE)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = alpha.1 ~ Treatment + Temperature, data = ulva)
## 
## $Treatment
##           diff        lwr         upr     p adj
## 1-0 -0.8506250 -1.1973752 -0.50387476 0.0000000
## 2-0 -1.1192708 -1.4660211 -0.77252059 0.0000000
## 3-0 -1.4839375 -1.8306877 -1.13718726 0.0000000
## 4-0 -1.5296250 -1.8763752 -1.18287476 0.0000000
## 2-1 -0.2686458 -0.6153961  0.07810441 0.2109076
## 3-1 -0.6333125 -0.9800627 -0.28656226 0.0000101
## 4-1 -0.6790000 -1.0257502 -0.33224976 0.0000018
## 3-2 -0.3646667 -0.7114169 -0.01791643 0.0338202
## 4-2 -0.4103542 -0.7571044 -0.06360393 0.0113476
## 4-3 -0.0456875 -0.3924377  0.30106274 0.9963020
TukeyHSD(ulva_photosyn_model_aov, "Temperature", ordered = FALSE)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = alpha.1 ~ Treatment + Temperature, data = ulva)
## 
## $Temperature
##            diff         lwr       upr     p adj
## 27-20  0.303925  0.07348064 0.5343694 0.0059250
## 30-20  0.252950  0.02250564 0.4833944 0.0275126
## 30-27 -0.050975 -0.28141936 0.1794694 0.8607823

alpha – Plot the results

plot(allEffects(run5_6_photosyn_model_noint))

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.